Flat Energy histogram version for Interacting Growth Walk 
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Interacting Growth Walks is a recently proposed stochastic model for studying the coil-globule 
transition of linear polymers. We propose a flat energy histogram version for Interacting Growth 
Walk. We demonstrate the algorithm on two dimensional square and triangular lattices by calcu- 
lating the density of energy states of Interacting Self Avoiding Walks. 
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Monte Carlo methods have emerged as a power- 
ful and reliable tool for simulating several complex phe- 
nomena in statistical physics, see, e.g., P, Conven- 
tional Monte Carlo methods [l| are found inefficient for 
simulating models with complicated energy landscapes, 
such as spin glasses or protein folding. There are many 
variants of the sampling techniques proposed to address 
these issues, see, e.g., parallel tempering Q. Several ef- 
forts have focussed on algorithms that aim to produce a 
nearly uniform distribution in one or more of the macro- 
scopic observables, such as energy or number of particles, 
within a predetermined range, which often fall in the class 
of so-called flat histogram methods, see below. 

Flat histogram techniques permit enhanced flex- 
ibility in sampling of energy space in that the system is 
given a greater probability for escaping low-lying energy 
minima as compared to traditional Boltzmann sampling. 
Multicanonical, transition matrix, and Wang- Landau al- 
gorithms [3| are some examples. These approaches have 
been used extensively to study liquid crystals, protein 
folding, and polymer phase behaviour |g] . This letter 
focusses on a flat energy histogram method to study the 
coil-globule transition of linear polymers. 

Long polymer chains in a good solvent have been 
studied extensively over the last five decades. Substan- 
tial progress has been made on elucidation of macroscopic 
properties as well as scaling behaviour of an isolated poly- 
mer chain from lattice models, see e.g. Mostly, these 
have been based on self avoiding walks, as they incorpo- 
rate in their basic definition excluded volume effect. 

Self Avoiding Walk (SAW) is a random walk in 
which a walk cannot visit a site more than once Q • In or- 
der to study the thermodynamic properties of SAW, one 
needs to assign energy to each conformation. A standard 
way of doing this for a lattice SAW is to assign an energy 
e to each nonbonded nearest neighbor (nbNN) contact in 



the walk. So, a conformation with m such contacts has 
energy E = m e. A SAW with energy assigned in this 
fashion is called an interacting self-avoiding walk(ISAW) 
*7:]. For the study of coil to globule transition of linear 
homopolymers, one can set e = — 1 without loss of gen- 
erality. 

Many conventional SAW Monte Carlo algo- 
rithms ^ and their variants have been formulated for 
obtaining equilibrium properties of ISAW. In a recent 
paper Rechnitzer and Rensberg Q define a random vari- 
able called 'atmosphere'. It is the number of possible 
ways a self avoiding walk of length n on a lattice can 
proceed to create a walk of length n + 1. They showed 
that from the statistics of atmosphere one can estimate 
the connective constant and entropic exponent [9]. 

More recently Prellberg and Krawczyk inter- 
preted the atmosphere as proportional to Rosenbluth and 
Rosenbluth (RR) weight [10| in a stochastic growth algo- 



rithm They show that RR walks in conjuction with 
the so-called 'dynamic' pruning and enrichment [l^ pro- 
vides a flat histogram Monte Carlo algorithm. 

In a recent work [l3l | , we have generalized the no- 
tion of atmosphere to other possible growth algorithms; 
we have shown that the average of atmosphere taken over 
ensembles of different stochastic growth algorithms give 
an estimate of the density of states of ISAW. In particular 
we have considered the Interacting Growth Walk (IGW) 
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When RR walks are kineticaly interpreted one 
gets Kinetic Growth Walk (KGW) [lB| ■ In the same way 



when the PERM-B walks of Grassberger [17| are kineti- 



cally intepreted one gets Interacting Growth Walk 

In IGW the first step from the origin to near- 
est neighbour site is taken randomly and independently 
with probability 1 /z where z is the coordination number 
of the lattice. Subsequent steps are taken with local IGW 
probability pj. Let us say that at j*'' step, Uj unvisited 
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nearest neighbour sites are available for the walk to pro- 
ceed. The probability for choosing the a*'' unvisited site 
is given by 



(1) 



where Pg is a model growth parameter; mnbNN(<^j) is 
the number of nonbonded nearest neighbour contacts the 
walk would make if it were to step onto the a^^ near- 
est neighbour site. If the number of unvisited nearest 
neighbour sites is zero at any stage, the walk is termi- 
nated (trapping); the entire walk is then discared and 
one starts all over again. For a given (3gi the probability 
of generating an N step IGW is given by, 



N 



pn = b{pg) X n 



(2) 



where B{I3g) is the normalization constant. The atmo- 
sphere for this N step IGW is given by (isj . 



^ 1 



(3) 



For = 0, every unvisited nearest neighbour 
site is chosen with equal probability. Therefore IGW re- 
duces to KGW in the limit /3g ^ [H]. For large 
/3g > the walks generated by the IGW are extremely 
compact. If /3g < 0, the walk would prefer to step onto 
that site that leads to lesser number of contacts. Hence 
for large negative /3g the walk generated by the IGW 
would be mostly extended as compared to those with 
/3g = 0. 

In this letter we propose a flat histogram Inter- 
acting Growth Walk. We show the atmosphere, when av- 
eraged over the flat histogram IGW ensemble, provides 
a powerful technique for estimating the density of states 
(DOS). 

In flat histogram IGW, the growth parameter Pg 
is taken as a random variable. The fluctuations of (3g are 
so adjusted that it ensures uniform accumulation of N- 
step ISAWs in all the energy bins. In other words if H{E) 
denotes the energy histogram of iV-step ISAWs, then the 
random variation of (Hq from one step to the next in a 
walk and from one walk to the other in the ensemble, 
leads to a flat H{E). 

The energy of an n step ISAW having m contacts 
is given by En — me. Let m™"^ denote the maximum 
number of contacts possible in an n step ISAW. Then the 
number of energy levels is m™"^ -|- 1. These are indexed 
by k = 1,2, ...(to™°^ + !)• Thus E„ can take values from 
zero to m^^^e and the correspoding energy histogram is 
represented by Hn{k). We also define an array Gn{k), 
called atmosphere sum, described below. 

In an A^-step flat histogram IGW algorithm, 
Hn{k) and Qn{k) are initially set zero for all n{l < n < 



N) and V A;. In a simulation run, if a walk at n*^ step has 
an energy E which corresponds to index fc, then Hn{k) 
and Gn{k) incremented as. 



where 



Hn{k) = Hn{k) + 1 

Gn{k) = Qn{k)+A^, 



1 

A„ = n- 



(4) 



(5) 



is the atmosphere [13[ for n*'*step for that run and pj 
is the probability for j*^ step which depends upon the 
value of Pg at that step. The implementation of the 
above algorithm proceeds as follows. 

The first step (n = 1) of the walk starts from the 
origin with probability l/z, where z is the coordination 
number of the chosen lattice. For n = 1, m™"^ — 
and Ai = z\ the energy histogram and atmosphere sum 
corresponding to the first step is incremented as -ffi(l) = 
Hi{l) -I- 1 and Qi{l) — Gi{l) + z. All subsequent steps 
(j > 2) are taken with probability pj given by Eq. ([T]). 
The value of /3g for each n(> 1) is chosen depending 
upon the difference between _ff„(?7-f 1) and Hn{ri), where 
T] = 1,2, ...(m™^ -I- 1). In order to span all energy levels 
the following procedure is adopted in this algorithm. 

Let 



A = 

S = 



min { I A 



Hniv)} e 



I , /3^"" } 



(6) 
(7) 



where ^ is a random number that takes values uniformly 
between zero to one and /3q"^ is the maximum value 
assigned to avoid overflow/ underflow problem. The value 
of /3g is chosen as 



f3G 



if A < 
if A > 
if A^O. 



(8) 



During each growth ste,p Pg takes a value between 
_pmax g^j^j^ (3^°-^ . The corresponding _ff„ and Qn are 
incremented as in Eq. ([?]). The walk is continued till it 
gets trapped or n reaches N . The above procedure is 
repeated for several number of Monte Carlo attempts M. 
Estimator for density of states, g{E), corresponding to 
the energy levels fc = 1, 2, ..(m™''^ -f 1) of ISAW V n < A 
is then given by 



9T{k) 



Qn{k) 
M ■ 



(9) 



From the density of states we can calculate the 
canonical partition function Z — (^(iS) exp(— /3i?) 
where (3 = l/fc^T. From the partition function one can 
calculate any thermodynamic property. 
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FIG. 1: Flat histogram IGW with Pa"'' = 10: Energy dis- 
tribution for A*' = 199 and 299. The energy distributions are 
reasonably flat. 



A general property of IGW is that with increase 
of Pg longer walks can be generated with increasing prob- 
abilities [l5|. In otherwords, attrition is less when /3g is 
large. In fact on a two dimensional square lattice, there 
is no attrition when /3g ~> 
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flat histogram IGW the parameter that control fluctua- 
tion of Pg is Pg""^ ■ Hence if Pq'^^ is large, flat histogram 
IGW can generate longer walks spanning uniformly all 
the energy levels. 

The flat histogram IGW described above was im- 
plemented on a two dimensional square lattice. Figure [1] 
depicts the energy histogram for N = 199 and 299 for flat 
histogram IGW with = 10. We find that proper 

choice of /3g would lead better statistics of density of 
states. Thus in our simulation on a square lattice we 
have chosen P^''"" = 1.0 and N < 300. The reduced his- 
togram hn{k) = Hn{k) / max{Hn{k)} Vfc for n = 25 to 
300 (insteps of 25) is shown in Figure [H The histograms 
are flat. The estimated density of states for shorter walk 
length compared with exact enumeration results is shown 
in Figure [H The simulation results match well with the 
exact results. The statistical error of each Monte Carlo 
data point does not exceed ±1%. The plot of estimated 
density of states versus energy for walk length n = 10 to 
120 (insteps of 10) is shown in Figured! These results 
were obtained over 10^ Monte Carlo attempts. 

We have calculated fluctuations of energy, 
a^{E) = {E'^) — (E)"^ as a function of p. This is shown in 
FigureH The value of /3 - (1/1.54) - 0.649 [7] at which 
phase transition is expected to occur is marked by verti- 
cal line. The fluctuations are maximum at the transition 
temperature. 

The "art" of making flat histogram IGW to es- 
timate exactly the DOS of ISAW for various lattices and 
various dimensions are based to a large extent on a suit- 
able choice of Pq°-^ . Flat histogram IGW simulation on 
triangular lattice also were carried out with P^"-^ = 0.8. 




FIG. 2: Plot of reduced histogram hn{En) versus energy En 
and length n = 25 to 300 (insteps of 25: right to left) for 
ISAW. These have been obtained from flat histogram IGW 
on a square lattice with = 1.0 and M = 10* attempts. 
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FIG. 3: Semilog plot of estimated density of states g'^^(En) 
versus energy En and length n = 5 to 25 (insteps of 5) for 
ISAW on a square lattices obtained from flat histogram IGW 
with = 1.0; M = 10*. The Monte Carlo resuhs (★) 

matches very well with exact (-I-) density of states. ( Monte 
Carlo errors are very small; maximum of error does not exeed 
one percent and are within the symbol size on the scale). 



The plot of reduced histogram and density of state are 
shown in Figure [5] . 

In conclusion we have presented a flat energy his- 
togram method to study Interacting Self Avoiding Walks. 
This method is based on IGW algorithm. Monte Carlo 
simulation of flat histogram IGW for longer self avoiding 
walks with large Pq°'^ one can obtain resonably flat en- 
ergy histogram. By optimizing Pq°-^ one can estimate 
density of states with good statics. We have carried 
out flat histogram IGW simulation of self avoiding walks 
on square and triangular lattices and presents results on 
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FIG. 4: Semilog plot of estimated density of states g'^^{En) 
versus energy En and length n = 10 to 120 (insteps of 10: 
right to left) for ISAW on a square lattice obtained from flat 
histogram IGW with /J^""^ = 1.0; M = 10^ 




FIG. 5: Energy fluctuation per walk length a^{E)/N versus 
inverse temperature /3 for ISAW on a square lattice for N = 
30(+),50(x),70(*) and 90(-). The vertical lines denote the 
expected transition temperature in the thermodynamic limit 
N ^ oo. 



DOS of ISAW and flat energy histogram. 
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FIG. 6: Flat histogram IGW on a triangular lattice, (a) Plot 
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various n; = 0.8; M = 10*. 
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